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Abstract 

We investigate the thermodynamic properties of a trapped Bose 
gas of Rb atoms interacting through a repulsive potential at low but 
finite temperature ( ksT < /u < Tc ) by Quantum Monte Carlo method 
based upon the generalization of Feynman-Kac method[l] applicable 
to many body systems at T=0 to finite temperatures. In this letter 
we report temperature variation of condensation fraction, chemical 
potential, density profile, total energy of the system, release energy, 
frequency shifts and moment of inertia within the realistic potential 
model( Morse type ) for the first time by diffusion Monte Carlo tech- 
nique. The most remarkable success was in achieving the same trend 
in the temperature variation of frequency shifts as was observed in 
JILA[2] for both m = 2 and m = modes. For other things we agree 
with the work of Giorgini et al[3], Pitaevskii et al[4] and Krauth[5]. 
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1 Introduction 



After the experimental realization [6] of Bose Einstein Condensation and mea- 
surement [7] of thermodynamic quantities in alkali gases, there has been a 
growing interest in the theoretical study of these system since these could 
easily be modelled as systems with weakly interacting condensate and along 
with other things theoretical predictions about thermodynamical properties 
were possible in terms of simple quantum statistical mechanics. Most of these 
studies involve computational techniques to solve the relevant many body sys- 
tem and based on mean field theory such as the Gross Pitaevskii[8] technique. 
Despite their success in explaining the ground state properties, predictions 
in finite temperature properties become only approximate and often can lead 
to incorrect predictions as we have seen from the discripencies[9-ll] between 
theory and experiment in explaining JILA top data . As a matter of fact 
mean field theory breaks down near Tc- It is therefore necessary to develop al- 
ternative computational methods which can solve these many body problems 
more accurately and rigorously. 

Thermodynamics of Bose gases was studied before at a higher 
temperature ( ksT » Tiu ) hy a, semiclassical treatment [3]. Since effects 
of interactions become more pronounced at low temperatures we restrict our 
discussions at low but finite temperature ( UbT < /i < Tc). At low temper- 
ature the de Broglie wavelength of the atoms become appreciable, the study 
of thermodynamic behaviour at low temperatures ( of the order of harmonic 
oscillator temperature ) requires a quantum description of a lowlying elemen- 
tary modes. As Quantum Monte Carlo technque and many body theory are 
closely connected, in this letter we present a quantum monte Carlo method 
namely Generalized Feynman-Kac method (GFK) [12,13] to study the ther- 
modynamic properties of a Bose gas. From the equivalence of the imaginary 
time propagator and temperature dependent density matrix, finite tempera- 
ture results can be obtained from the same zero temperature code by running 
it for finite time. The first Monte carlo calculations [5] on BEG deals with 
temperature dependence of condensation fraction and the other remarkable 
Monte Carlo calculation[14] deals with the ground state properties. We cal- 
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culate temperature variation of condensation fraction, total energy, release 
energy, frequency shift, chemical potential and moment of inertia for system 
of 100 RlP''^ atoms. 



2 Theory 

All the thermodynamic quantities of interest are connected with the evalua- 
tion of eigensolution and eigenencrgy of the the many body system. In this 
non mean field approach, we consider the full Hamiltonian for 100 R}/'^ atoms 
interacting through Morse potential at low but finite temperature and the 
solution of corresponding Schroedinger equation in path integral representa- 
tion. We first describe the T=0 version of the GFK formalism[12,13] and then 
generahze it to finite temperatures. For the Hamiltonian H — —A/2 -|- V{x) 
consider the initial value problem 

= {-- + v)uM 

u{0,x) = f{x) (1) 

with X E and u{0,x) — 1. The solution of the above equation can be 
written in Feynman-Kac representation as 

u(t, x) = E^exp{- r V(X(s))ds} (2) 
Jo 

where X(t) is a Brownian motion trajectory and E is the average value of the 
exponential term with respect to these trajectories. To speed up the conver- 
gence one can incorporate importance sampling in the underlying stochastic 
process and the lowest energy formula for eigenvalue for a given symmetry 
obtained from the large deviation priniciple of Donsker and Varadhan [15] 
can be written as 

1 /■* 
A = At - lim -lnE^exp{- / Vp{Y{s))ds} (3) 

>oo t Jo 

where Y(t) is the diffusion process which solves the stochastic differential 
equation 
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The temperature dependence comes from the reaUzation that the 
imaginary time propapagator is identical to the temperature dependent den- 
sity matrix \lt ^ (3 — 1/T holds. 

This becomes obvious when we consider the eqs[16] 



dt2 

and 



H2k{2, 1) (5) 



-^^H2P{2,1) (6) 

and compare 

fc(2,l)=^0,(a;2)0.*(a;i)e-(*-*^)^^ (7) 

i 

and 

p{2,l)^J2<l>,{x2W{x,)e-^'=^ (8) 

i 

For Zero temp FK we had to extrapolate to ^ =^ oo. For finite run time t 
in the simulation, we have finite temperature results. Here we show how to 
change our formalism to go from zero to finite temperature. We begin with 
the definition of finite temperature. A particular temperature 'T' is said to 
be finite if AE < kT holds. The temperature dependent density matrix can 
be written in the following form 

p{x,x',l3) = p('^)(x,x',/3)x < exp[- fVp[X{s)]ds\ >drw (9) 

At finite temperature thus free energy can be written as 

F = -lnZ{x,p)/p = -lnZ\x,l3)/l3-ln < exp[- [\p[X{s)]ds] >drw/P 

J 

(10) 
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2.1 The stationary state of the condensate 

To calculate the condensate energy and condensate density, we consider a 
cloud of N atoms interacting through repulsive potential and placed in a 
three dimensional anharmonic oscillator potential. We will assume that the 
condensate system has relaxed to a stationary state and at low energy the 
time independent Schroedinger equation representing the stationary state 
can be written as 

1 ^ 

[-A/2 + Vint + ^ + Vi^ + = Eil^^i^ (11) 

^ 1=1 

where \ Y,f^i[xi^ + + \z^] is the anisotropic potential with anisotropy 
factor A = ^. Now 

Vint = VMorse = E^inj) = ^ D[e-"('-'-«) (e-^^'-''^) - 2)] (12) 

Here we assume that the condensate oscillates in a static thermal bath. There 
is no interaction between the condensate and the thermal bath. The principal 
effect of finite temperature on the excitations is the depletion of condensate 
atoms. In the dilute limit and at very low energy only binary collisions 
are possible and no three body recombination is allowed. In such two body 
scattering at low energy first order Born approximation is applicable and the 
interaction strength 'D' in the dimnsionless form(7) turns out to be 

3 

_ 4 ^ ^ ^ 4 g ^ /^g) 

_ggaro(-garo _ 16) ^ ^ 

For more details, one should look at [17]. In the above expression a is the 
scattering length of Rb, a is the width of the Morse potential, Tq is the 
minimum of the potential well, 's' is the length scale in units of harmonic 
oscillator and A is the anisotropy factor. Here we have used [17] a = .29 
in harmonic oscillator units, tq = 9.67 in harmonic oscillator units, a = 
52 X 10-1° m, s = .12 x 10"^ m and A = \/8 For more details, one should 
look at [17]. 
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2.2 Effect of noncondensate 

In the case of noncondensate the system can be considered as a thermal gas. 

To calculate noncondensate energy and density we need to study the effect of 
noncondensate exphcitly and consider the following stationary state for the 
thermal gas. 

[-A/2 + 2Vir,t + Vtrap]lpj{T) = E^,iP^{t) (14) 
1 ^ 

[-A/2 + 2Vint + - E[^^' + + A;2,2]^,-](f) = ^„e^,(f) (15) 
^ 1=1 

The basis wavefunction which describes the noncondensate should be 
chosen in such a way that it is orthogonal to ■00 ^ in Eq.(ll) The most 
common way to achieve a orthogonal basis in Schroedinger prescription is 
to consider the dynamics of noncondensate in an effective potcntial[18,19] 
Veff = Vtrap + 2Vi„(. Thc factor 2 represents the exchange term between two 
atoms in two different states. The energy in the case of lowest lying modes 
then corresponds to E — Ec + E^c- One can calculate the E^c using the same 
parameters as discussed in Sec 2.1 
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3 Results and discussions 



3.1 The effect of temperature on frequency shifts 
3.1.1 ixi—2 mode 

Using Eq.(3), one can calculate the lowest lying energies due to any symme- 
try. To calculate the frequencies for m = and m = 2 modes, one needs 
to find the energy differences of each of these states and the ground state. 
Underneath we show the data for frequency shifts for both m = 2 and m = 
modes. For m = 2 mode, considering motion of condensate only we achieve 
the downward shift of data all the way to T = 0.9Tc. But for the m = 
mode we need to consider the dynamics of thermal cloud also as discussed 
above. 
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Figure 1: Effects of temperature on m=2 mode; this work. The top cuve from 
equivalent T=0 system[method 2], the bottom curve by putting temperature 
directly [method 1]. Both show agreement with JILA experimental data[2] 
all the way up to 0.9Tc and Ref[20,21] 
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3.1.2 m=0 mode 

For m = mode, considering the dynamics of condensate ( Sec 2.1) alone 

we do not get the upward shift as observed in JILA experimental data. But 
when the motion of thermal cloud (Sec 2.2) is considered in a dynamical 
manner, we observe the expected upward shift at around T = 0.7Tc. 




Figure 2: Effects of temperature on m=0 mode from GFK considering 
noncondensate dynamics [this work], shows resemblence with JILA [2] and 
Ref[20,21] 

3.2 The effect of temperature on condensate density 

The condensate density can be evaluated solving Eq.(ll) and using Eq.(9). 
Underneath we plot the axial density due to condensate along x axis. 
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Figure 3: Axial density profile due to Condensate at temperature T=0.48 




Figure 4: Axial density profile due to condensate at temperature T=0.6 

Fig 3, 4 and 5 are concerned with condensate density profile vs x 
for temperatures T — 0.48, T — 0.6 and T — 0.7 respectively for a system 
of 100 Rb atoms. We see that the center density for condensate increases as 
temperature is increased. This agrees with the earlier work of Krauth[5] 
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Figure 5: Axial density profile due to condensate at temperature T=0.7 
3.3 The effect of temperature on the total density 



Following the theory in Sec 2.1 and 2.2 and using Eq.(9) one can calculate the 
total density of the Bose gas. Underneath we show the temperature variation 
of the total density along the x axis. 




X 

Figure 6: total density profile at temperature T=0.48 
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Figure 7: total density profile at temperature T=0.6 



0.000278 




Figure 8: total density profile at temperature T=0.7 

On the other hand in Fig 6, 7 and 8 which represent the total 
density profile vs x at T = .48,7" = 0.6 and T — 0.7 respectively, we see the 
opposite trend as the center density decreases as the temperature is raised. 
This also is in agreement with the ealier observation by Krauth[5] 
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3.4 Effect of temperature on the total energy of the 
Bose gas 

Total energy of the system The condensate and noncondensate energies 
can be evaluated solving Eq. (11) and Eq.(15) and using Eq(3). Then com- 
bining condensate and noncondensate energy we get the total energy of the 
system. The total energy of condensed and noncondensed component of a 
trapped Bose gas is a combination of E^o and Ei^t for each component 
separately. The trend in the temperature variation of total energy is found 
to be the same as in [3,4]. In principle, specific heat can be calculated as the 
temperature derivative of total energy/particle keeping the confining poten- 
tial constant. We will report it for bigger systems somewhere else. Release 
energy can be represented defined to be the energy obtained after switching 
off trap. We get similar trend the temperature variation of release energy as 
we see in ref [4] 
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Figure 9: Total energy /particle as a function of reduced temperature. Inset: 
Release energy as a function of temperature 
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3.5 Effects of temperature on condensation fraction 

Interaction lowers the condensation fraction (-^ ) for repulsive potentials. 
Some particles always leave the trap because of the repulsive nature of the 
potential and on the top of it, if temperature is increased further, more par- 
ticles will fall out of the trap and get thermally distributed. This decrease in 
condensation fraction eventually would cause the shifts in the critical tem- 
perature ( Tc decreases ). We would observe this in (Fig. 10 ). Earlier this was 
done by W. Krauth[5] for a large number of atoms by path integral Monte 
Carlo methc^ 




Reduced Temperature T/T^ 

Figure 10: Condensation fraction vs Reduced Temperature ; this work. The 
innermost curve corresponds to the 100 interacting atoms and the outermost 
curve corresponds to the noninteracting case. The middle curve corresponds 
to the interactive case in ref[3]. The number of condensed particles decreases 
with the interaction 
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3.6 Effects of temperature on Chemical potential 

The chemical potential /j. can be written in terms of different contributions 

to the energy, namely -Efem, Eho and Eint as follows [4]. 
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Figure 11: Chemical Potential vs Reduced Temperature 
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3.7 Effects of temperature on moment of inertia 

Moment of inertia The deviation of moment of inertia from its rigid value 
is given by the useful expression [22] 



e 

9~~N'o<r^>o +Nt < r2 >t 



(17) 



where <>o and <>t denote the average taken over the condensate and 
noncondensate densities of the Bose gas respectively. 




Figure 12: Moment of inertia 9 divided by its rigid value 9rig as a function 
of T/T, 



©rig 



At T = 0, Nt = and ^ = 0. On the other hand f o T = T^, 
1 as A^o = 0. 
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4 Conclusions: 

For the first time we have calculated finite temperature properties beyond 
meanfield approximation( GP etc) by Quantum Monte Carlo technique. We 
have calculated spectrum of Rb gas by considering realistic potentials like 
Morse potential etc. instead of conventional pseudopotentials for the first 
time. Since a dilute gas consisting of atoms is a bosonic system the 
random walk in GFK is exact in the limit scale, time for walk and the num- 
ber of walks get arbitraly large and it turns out that the this method is a 
potentially good candidate for a sampling procedure for Bose gases at all 
temperatures. 

We are dealing with only 100 atoms. Nonetheless we have been 
able to show the variation of condensate and noncondensate density with 
temperature as a hallmark of BEG and lowering of condensation fraction in 
the case of interacting case which is very unique with a system of trapped 
atoms compared to uniform system. 

In our non mean field study for finite temperature excitaion spec- 
trum for m=2 mode, we see agreement with experimental study all the way 
to T = 0.9Tc [Fig. 1]. For m = mode, considering the motion of thermal 
cloud explicitly, we observe the upward shift of data[Fig 2] in JILA[2] and 
Ref[20,21]. Since at low temperature we solve the many body problem by 
considering the full Hamiltonian with realistic potential full quantum me- 
chanically and nonperturbatively, obviously the above modes are collective 
in nature and correspond to the m=2 and m=0 experimental modes[7]. 

Since we are dealing with very small number of particles we cannot 
compare our data with other existing results directly on the same graph. 
We can neither use scaling method since we are working on a small system 
at a low temperature regime ks < jJ. < T^. We can study and compare 
the temperature variation of different quantities such as total energy, release 
energy, moment of inertia with those in the literature [3, 4, 5] only qualitatively. 

The method is extremely easy to implement and our fortran code 
at this point consists of about 270 lines. We employ an algorithm which is 
essentially parallel in nature so that eventually we can parallelize our code 
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and calculate thermodynamic properties of bigger systems with the scaling 
property ( of the order of 2000 atoms ) taking advantage of the new computer 
architectures This work is in progress. We are continuing on this problem 
and hope that this technique will inspire others to do similar calculations. 
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